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C^ • Abstract We present a discretization for Darcy's problem using the recently devel- 

oped Mimetic Spectral Element Method [19]. The gist lies in the exact discrete rep- 
resentation of integral relations. In this paper, an anisotropic flow through a porous 
medium is considered and a discretization of a full permeability tensor is presented. 
The performance of the method is evaluated on standard test problems, converging 

r"| ' at the same rate as the best possible approximation. 
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1 DARCY FLOW 



Anisotropic heterogeneous diffiision problems are ubiquitous across different scien- 
tific fields, such as, hydrogeology, oil reservoir simulation, plasma physics, biology, 
etc [10]. Darcy's equation describes a steady pressure-driven flow through a porous 
t~^ , medium where fluxes and pressure are linearly related. 
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u - grad p-OmQ (la) 
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q = Ku in^ (Ic) 

q = qo indQ (Id) 



j^ \ where u is the fluid velocity, p the pressure, q the mass flux and (p the prescribed 



source term. Without loss of generality let the viscosity, ju = 1, and consider a per- 
meability symmetric, positive definite tensor denoted by K. 

In a three-dimensional setting are: four types of submanifolds {points, lines, sur- 
faces and volumes); and two orientations {outer and inner, as an example see Fig- 
ure 1). Tessellation divides the physical domain in a set of these geometric objects to 
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which we associate discrete variables, i.e. integral quantities. Thus, associated with 
every physical variable is a correspondent geometric object, this symbiotic relation 
between physics and geometry is the core of mimetic methods. Many scholars are 
aware of this relationship [3, 5, 7, 23]. 

Outer orieatarion Imier orieatarion 




Fig. 1 Consider a line wliere we can liave two types of orientation: Outer - around tlie line; Inner 
- along the line. 



Starting from the mass balance equation, (lb), 

I divqdy= ( qnd5= \ (pdV, (2) 

Jv JdV Jv 

it is clear that the divergence in a volume is equal to the sum of the surface integral 
quantities, i.e. oriented fluxes. Thus, we will associate mass fluxes, q, with quantities 
that go through surfaces. This equation therefore tells us that the right hand side term 
is associated to outer-oriented volumes. 

Similarly, using Newton-Leibniz relation for equation (la). 



I grad/?dC= | p^p(B)-p(A)^ (udC, 

Jc JdC Jc 



(3) 



the fluid velocity, u, is represented along lines and p is represented by the values in 
points. From (3) we deduces that u and p are inner-oriented variables. 
The constitutive/material relation relation (Ic) is given by, 

q = Ku, (4) 

which defines how quantities associated to inner-oriented lines relate to quantities 
associated to outer-oriented surfaces. Whereas equation (2) and (3) can be exactly 
satisfied on a finite grid the constitutive equation (4) needs to be approximated. 

The importance of respecting the geometric nature in physics is discussed in [9]. 
Figure 2 summarizes the geometric character of the Darcy's problem. 

We will denote the space of variables associated to outer-oriented ^-dimensional 
objects by A^ (A1) and the space of variables associated to inner-oriented fe-dimensional 
objects by A^ (M) as indicated in Figure 2. 

In this paper we will make use of the spectral element method described in [9, 1 9] , 
application of these ideas to Stokes' flow see [16-18]; Poisson equation for vol- 
ume forms [22]; advection equation [21]; derivation of a momentum conservation 
scheme [24]. Extension to compatible isogeometric methods see [11,12]. For appli- 
cations of these ideas in a finite diff'erence setting see Brezzi et al. [4]. In the context 
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Fig. 2 Darcy's flow problem geometric characterization. Fluxes, q, are associated with outer ori- 
ented surfaces; ^, is associated with outer oriented volumes; velocity, u, is associated with inner 
oriented lines; pressure, p, is associated with inner oriented points 



of finite element methods Arnold, FaUc and Winther [1] proposed a Finite Element 
Exterior Calculus. In a more geometric spirit Desbrun et al. [6] and Hirani [13] de- 
veloped the discrete exterior calculus (DEC). An application of the latter to Darcy 
flow can be found in [14]. 



2 DISCRETIZATION OF EQUATIONS 

In this section we will describe the discretization by defining the weak formulation. 
The approached followed here is similar to [2, 15]. 



For vectors associated with outer-oriented surfaces, A^(A1), we define the 
weighted inner product 

(a,b)AiK= r aK-'bdV. (5) 

Furthermore, we define bilinear maps ((■, 0)^ : A' (A1) x A"~^ (Al) — » R 

((u,v))ai:= r u-vdV. (6) 

and ((■, •))ai,k : A*^ (M) x A""*^ (At) -» R given by 



((u,v))yM,K 



JM 



'vdV. 



(7) 



For q e A^ (Al) and /? e A" (Al) and homogeneous boundary values we have 

((divq, p))f^ = -((q,gradp));^ 

= -((q,K-'[Kgrad/,]))^ (8) 

= ((q,gradip))^_^. 
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It is possible to define a new gradient operator, 

grad^ = -Kgrad. (9) 

Mixed formulation 

Starting from (1) and making use of the bilinear maps defined above we have for all 
vectors r E A"~' (A1) associated to outer-oriented surfaces 

((T,u-grad/?))y^, = 
((T,Ku-Kgrad/,))^_j, = 

(9) 

^^ (10) 

((t,Ku))ai,k + (r.grad^p)^ j^ = 

(ic) and (8) 

('r,q)M,K + ((divT, /?));k=0 

The constitutive equation is included in the last step by converting the bilinear form 
to a weighted inner product on A^ (At) as defined in (5). For ( lb) we take the bilinear 
map for the divergence of a vector q associated with outer-oriented surfaces and an 
arbitrary scalar function defined in inner-oriented points, y e A** (At), 

((divq,r))Ai = ((0'r))M- (11) 



The mixed formulation becomes: Find (q,/?) 


e[A^{M)y 


■ A^(M)}, 


given e 


A^ (M), for all (r 


,r) e [a^ (M) X a" (M)} such that. 








(T,q)Ai,K + ((divT, p))m-- 


--Q 




(12) 




((divq,y))^ -- 


--((<f>,7))M- 




(13) 



2.1 Basis functions 

For the high order representation we use Lagrange, /, (0, and edge functions, e, (^). 
Lagrange polynomials interpolate nodal values. The edge functions, derived by Ger- 



Mixed Mimetic Spectral Element method applied to Darcy's problem 5 

ritsma [8] are constructed such that when integrating over a Hne segment it gives one 
for the corresponding element and zero for any other line segment, 






(14) 



The relation between the Lagrange and the edge functions is given by, 



i-i 



eiiO^eii^d^, with e,.(^) = -^^. 



k=0 



dlk 
'di' 



Note that this definition implies 



dU 
d^ 



= e,(^)-e,+i(^). 



(15) 



(16) 



Extension to the multidimensional is obtained by means of tensor products. For 
more details see [19]. 



2.2 Mimetic discretization in 2D 



Expansion of unknowns in M? 

Let qeA^ (M) be expanded as, 

and the pressure, ph e A" (M) as. 
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(17) 



(18) 



Discrete divergence in 9? 

The divergence of q/, is then given by 



N N 



divq, = XZfe-^^-iJ + ^L-^L-i)^'^^^^-'-^''^ 



(19) 
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where we repeatedly used (16). The scalar e A^ (At) associated with outer-oriented 
volumes is expanded as 



N N 

0/1 = ^ ^ 4>i.jei (0 ej (77) . 
Equating (19) and (20) yields 

N N N N 

,=1 ;=1 ,=1 j=l 

[0] = E(2.i)[q]. 



(20) 



(21) 



(22) 



We see that the basis functions cancel from this relation. The matrix E*^^''^ relates 
the fluxes q'f . and q^. . to the volume integral 0,,,, as depicted in Figure 3. This fully 

discrete equation is a restatement of the integral relation (2). The matrix E^^''^ only 
contains the values -1,0 and 1 and is fully determined by the grid, see [19]. This is 
an incidence matrix showing the topological nature of the discrete divergence. 
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Fig. 3 Discrete representation of the action of the divergence in W? and I 



If we insert the expansions of our unknowns in (12) and (13) we obtain in R" 
the saddle point problem given by. 
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^^«-i) (e("-"-1))^mW 



j\^(n)£(n,n-l) 











(23) 



where M^''^ is the symmetric mass matrix obtain from the biUnear pairing 
between variables associated with outer-orientation and inner-orientation, (6), 



^(n-l) 



Mjj is the mass matrix obtained from the weighted inner product (5) and 
^(m.m-I) jjjg incidence matrix which relates fluxes over surfaces to volumes. 
The resulting system (23) is symmetric. 



The pressure which is represented on an inner-oriented grid (which is not explic- 
itly constructed in this single grid approach) is pre-multiplied by M^"' to represent 
it on the outer-oriented grid. 



3 NUMERICAL RESULTS 

The method derived in this paper respects the geometric nature of the problem. 
However, it is crucial to verify the numerical benefits of this approach. This section 
presents /ip -convergence studies for anisotropic permeability. 



3.1 Manufactured solution - Anisotropic permeability 



The first test case assesses the convergence for h- and /^-refinement of the mixed 
mimetic spectral element method applied to the Darcy model. This is a benchmark 
problem presented in [15]. The problem is defined on a unit square, Q - [-1, 1]^, 
with Cartesian coordinates with permeability given by. 



2 1 
1 2 



and the right hand side, (f> e L^ (M) given by, 

cf>^^^ = 2 ( 1 + x^ -H xy + y^) e^' dxdy. 
This results in an exact solution for pressure /? e A" (M) given by. 



(24) 



(25) 



(26) 



Figure 4 shows the h- and /? -convergence for the pressure in straight mesh. For the 
/i-convergence the expected rate of convergence is of (p + l), where p is the poly- 
nomial degree. The solid line for the interpolation error in the /^-convergence plot 
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is the L^-error from interpolating the exact solution, the solution converges expo- 
nentially. Both the numerical solution and the interpolated exact solution converge 
exponentially. 
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Fig. 4 Plots of the h- and p-convergence for anisotropic permeability given in (24). 



3.2 Layered medium 

A classical benchmark for Darcy flow codes is the piecewise constant permeability 
in a square [20]. Such a medium is called layered medium. 



Oa 



ifi<3'<i 



0.3 
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0.5 if3'>f 



(27) 



The fluid comes into the domain from the left to the right. Since the pressure depends 
linearly on x, horizontal constant velocity is expected in each layer. Figure 5. 
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Q_|. First of all the authors would like to acknowledge the reviewers for their 

•^ ' feedback. 

<N ■ 1 Reviewer 1 :: IC012-2-12.Rl.pdf 



.^ ' • An extra figure was added for clarity of what is meant with two orienta- 

tions (Figure 1). The reader is also directed to the paper on this issue, The 
geometric basis of numerical methods, where a careful analysis of geometry 
is given. 



Formula (14) is corrected. 



• Formula (27) is corrected since pressure belongs to A'^ (A4) the term dxdy 
should vanish. 

• Reference 1 was completed together with some other references that were 



> 

■^ ■ recently published. 

t^^ i English was carefully reviewed and includes all the comments. 

^ : 2 Reviewer 2 :: IC012-2-12.R2.pdf 

• A reference to the research areas and the FVCA6 benchmark is added to 
the first sentence. 

X- 

^ , • Citations referring to the work of Bossavit and Tonti were added. 



• Citations about Finite Element Calculus, Arnold et al; Mimetic Finite 
Difference, Brezzi et al; Analysis of Compatible Discrete Operator Schemes 
for Elliptic Problems on Polyhedral Meshes, Bonelle et al.; DEC on Darcy, 
Hirani et al. 

• In our formulation we avoid the term scalar product because of entries 
in this bilinear form come from different spaces, inner and outer oriented 
variables associated with different geometric objects. This is the main 
idea in this paper, therefore we prefer to use bi-linear form. 

• The term adjoint was removed. 



• The sentence was changed according to the reviewer suggestion. 
All the minor reviews were taken into account. 



